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We investigate the dynamics of the two-dimensional periodic Kolmogorov flow of a viscoelastic 
fluid, described by the Oldroyd-B model, by means of direct numerical simulations. Above a critical 
Weissenberg number the flow displays a transition from stationary to randomly fluctuating states, 
via periodic ones. The increasing complexity of the flow in both time and space at progressively 
higher values of elasticity accompanies the establishment of mixing features. The peculiar dynam- 
ical behavior observed in the simulations is found to be related to the appearance of filamental 
propagating patterns, which develop even in the limit of very small inertial non-linearities, thanks 
to the feedback of elastic forces on the flow. 

PACS numbers: 47.27.E-, 47.57.Ng 



I. INTRODUCTION 

It is well known that the addition of small amounts of 
long chain polymers can have a strong impact on flowing 
fluids' properties. One of the most remarkable effects of 
highly viscous polymer solutions which has been recently 
observed in experiments is the development of an "elastic 
turbulence" regime in the limit of small Reynolds num- 
ber and strong elasticity [l| . The flow in this regime dis- 
plays a rich dynamics, of the type usually encountered 
in high Reynolds number turbulent flow in Newtonian 
fluids. Typical features of turbulent flows (e.g., broad 
range of active scales, growth of flow resistance) are ob- 
served in the experiments by Steinberg and collabora- 
tors Q, even at low velocity and high viscosity, i.e. 
in the limit of vanishing Reynolds number. As a con- 
sequence of turbulent-like motion, elastic turbulence has 
been proposed as an efficient technique for mixing in very 
low Reynolds flows, such as in microchannel flows [3|-[5j. 
Despite its great technological interest, elastic turbulence 
is still only partially understood from a theoretical point 
of view. Presently available theoretical predictions are 
based on simplified versions of viscoelastic models and on 
the analogy with magneto-hydrodynamics (MHD) equa- 
tions d, 0- In a previous work [H we showed, by means 
of numerical simulations, that the basic phenomenology 
of elastic turbulence is reproduced by a simple model of 
polymeric fluid and in an idealized geometry, namely a 
two-dimensional periodic shear flow without boundaries. 
Notwithstanding quantitative differences, most observed 
features have a strong qualitative resemblance with ex- 
perimental results. 

The mechanism through which polymer molecules can 
influence properties of a flow is their extreme extensibil- 
ity. Polymers, typically composed by a large number of 
monomers, at equilibrium are coiled in a ball of radius 
Ro. In presence of a non- homogeneous flow, the molecule 
is deformed in an elongated structure characterized by its 
end-to-end distance R which can be much larger than Rq. 



The deformation of molecules is the result of the competi- 
tion between the stretching induced by velocity gradients 
and the entropic relaxation of polymers to their equilib- 
rium configuration. Experiments with DNA molecules 
0, [Toj| show that this relaxation is linear, provided that 
the elongation is small compared with the maximal ex- 
tension R <C Rmax, and therefore can be characterized 
by a typical relaxation time r . Besides the Reynolds 
number i?e, commonly used in hydrodynamics, in the 
study of flowing polymer solutions, it is useful to intro- 
duce the Weissenberg number Wi. This is defined as the 
product of r and the characteristic deformation rate, and 
it measures the relative importance between polymer re- 
laxation and stretching. When Wi <C 1 relaxation is fast 
compared to the stretching time and polymers remain in 
the coiled state. For Wi ^> 1, on the contrary, polymers 
are substantially elongated. The transition point is called 
the coil-stretch transition and occurs at Wi = 0(1). For 
Weissenberg numbers much larger than unity, polymers 
can provide a significant feedback effect on the dynami- 
cal behavior of the fluid they are suspended into, so that 
they behave as active objects capable of modifying the 
properties of the flow. 

Because of the feedback effect, above the coil-stretch 
transition, a laminar flow can become unstable even in 
the small Re limit of vanishing inertial non-linearities. 
Purely elastic instabilities have been studied in different 
flow configurations. In a curvilinear shear flow as, e.g., 
Couette flow between rotating cylinders, the driving force 
of the elastic instability is the "hoop stress" (see, for in- 
stance, [12| page 62), a volume force in the direction of 
curvature due to the first normal stress difference. The 
mechanism for this instability was first proposed in [l3| 
and experimentally verified in [14]. In the framework of 
flows with curvilinear streamlines, a dimensionless cri- 
terion for the critical onset conditions, which applies to 
various flow geometries, was established in [HI, [l6|. The 
role of flow-microstructure coupling in pattern formation 
in the Taylor- Couette configuration was also examined in 
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numerical simulations of the 3D FENE-P (finitely exten- 
sible non-linear elastic-Peterlin) model [17]. On the other 
hand, purely elastic instabilities in a parallel flow with 
rectilinear streamlines were obtained for the first time in 
a 2D Kolmogorov flow of a viscoelastic fluid described by 
the linear Oldroyd-B model [18|- Also for a more realis- 
tic shear flow of a viscoelastic fluid in a straight channel, 
there exist theoretical results exhibiting the elastic in- 
stability [19j, although they are still not experimentally 
verified. Elastically induced pattern formation and pro- 
duction of fluid mixing were also studied in extensional 
flows, in experiments [20], and numerical simulations of 
a Stokesian linear viscoelastic model [21] . 

In this paper, we study the dynamics of dilute polymer 
solutions by means of direct numerical simulations of a 
general linear viscoelastic model, in a flow configuration 
corresponding to the two-dimensional Kolmogorov shear 
flow. Our simulations for values of the stability parame- 
ters (Re and Wi) close to the elastic instability threshold, 
exhibit a transition to new intriguing flow states, in the 
form of "elastic waves" , driving the establishment of mix- 
ing features. In particular, at growing the elasticity of the 
solution, after a first transition to steady uniformly trav- 
eling patterns with different symmetry with respect to 
the reference laminar fixed point, we detect the appear- 
ance of temporal oscillations and, eventually, the onset of 
a strongly non-linear multiscale state disordered in both 
space and time. 

The article is organized as follows: in the next sec- 
tion we introduce the viscoelastic model. In section UlI Al 
we present numerical results concerning the transition to 
mixing states in terms of an analysis of the behavior of 
global observables. In section IIIIBI we show the results 
of numerical measurements concerning the spatiotempo- 
ral structure of the patterns emerging above the onset of 
the instability. Finally, section IIVI summarizes our con- 
clusions. 



II. MODEL 

To describe the dynamics of a dilute polymer solution 
we adopt the well known linear Oldroyd-B model [12| 

d t u + (u • V)u = -Vp + vAu + — V • cr + /, (1) 

r 

(),a + (u ■ V)a = (Vu) T cr + cr- (Vu) -2 ^—^ , (2) 

r 

where u is the incompressible velocity field, cr the sym- 
metric positive definite matrix representing the normal- 
ized conformation tensor of polymer molecules and 1 is 
the unit tensor. The trace trcr is a measure of local poly- 
mer elongation. The solvent viscosity is denoted by v and 
r] is the zero-shear contribution of polymers to the total 
solution viscosity v t — 1/(1 + 77) an d is proportional to the 
polymer concentration. In absence of flow, u = 0, poly- 
mers relax to the equilibrium configuration and cr = 1. 



The extra stress term V • a in (pQ) takes into account 
the elastic forces of polymers and provides a feedback 
mechanism on the flow. 

In order to concentrate on the intrinsic dynamics of 
the viscoelastic solution we use the simple geometrical 
setup of the periodic Kolmogorov flow in two dimensions 
|22| . With the forcing / = (Fcos(y/L), 0), the system 
of equations (pj2]) has a laminar Kolmogorov fixed point 
given by 

^0 = (U cos(y/L),Q) 

(3) 

_ / 1 + r 2 S sm 2 (y/L) -r Sf. sin (y/L) \ 
V -r^ sin (y/L) 1 / 

with F = [vU (l + r]))/L 2 [18]. The laminar flow fixes a 
characteristic scale L, velocity Uo and time T = L/Uq. 
In terms of these variables, one can define the Reynolds 
number as Re = and the Weissenberg number as 

Wi = 1 j R . The ratio of these numbers defines the elastic- 
ity of the flow El = Wi/ Re. Notice these are a-priori es- 
timates of the non-dimensional parameters; in the follow- 
ing we will refer to them as Reo and Wio- In the presence 
of velocity fluctuations superimposed to a non-zero mean 
flow [/, it is also possible to compute a-posteriori non- 
dimensional parameters as Re = ^ and Wi = ^ , once 
the average velocity U has been computed; of course, in 
general U 7^ Uq. 

It is well known that the Kolmogorov flow displays 
instability with respect to large-scale perturbations, i.e. 
with wavelength much larger than L. In the Newtonian 
case, the instability arises at Re c = V2 [23]. The pres- 
ence of polym ers changes the stability diagram of laminar 
flows 124, 251 and can induce elastic instabilities which are 
not present in the Newtonian limit [l3|, [3, EH, , even 
in the case of periodic flows (27j . In this respect, the 
Kolmogorov flow is no exception, and recent analytical 
and numerical investigations have offered the complete 
instability diagram in the (Wi, Re) plane [l8j . For the 
purpose of the present work, we just have to recall that 
linear stability analysis shows that for sufficiently large 
values of elasticity, the Kolmogorov flow displays purely 
elastic instabilities, even at vanishing Reynolds number 
(see Fig. 1 of [18]). Note that in this case the direction of 
the most unstable mode is at a small angle with respect 
to the basic flow, at variance with the purely hydrody- 
namic transverse instabilities. In the Newtonian case, 
non-transverse instabilities are observed for less symmet- 
ric basic flows only [28|- We also remark that the fact 
that the basic flow has rectilinear streamlines does not 
exclude the onset of the elastic instability [l9| . Above 
the elastic instability the flow can develop a disordered 
secondary flow which persists in the limit of vanishing 
Reynolds number and eventually leads to the elastic tur- 
bulence regime [29|. 
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III. ANALYSIS AND RESULTS 

The equations of motion (P)-© are integrated by 
means of a pseudo-spectral method implemented on a 
two-dimensional grid of size Lq = 2tt with periodic 
boundary conditions at resolution 512 2 . It is well known 
that the integration of viscoelastic models is limited by 
numerical instabilities associated with the loss of posi- 
tiveness of the conformation tensor [30| . These instabili- 
ties are particularly important at high elasticity and limit 
the possibility to investigate the elastic turbulent regime 
by direct implementation of equations (PQ)-(|2]). To over- 
come this problem, we have implemented an algorithm 
based on a Cholesky decomposition of the conformation 
matrix that ensures symmetry and positive definiteness 
|31[. While this allows us to safely reach high elasticity 
it has a price to pay in terms of limited resolution and 
large computing time. In the following, we will present 
the numerical results obtained using this approach. 



A. Transition to mixing states 

In order to investigate the transition to mixing states, 
we chose to work in a region of parameters close to the 
onset of the purely elastic instability. Here we fix Re = 1 
and increase the Weissenberg number from Wi = 6 to 
Wi > 20. The initial condition is obtained by adding a 
small random perturbation to the fixed point solution (j3j) 
and the system is evolved in time until a (statistically) 
steady state is reached. 

We start the analysis by considering the behavior of 

= /L,|2\ 



global quantities, such as the kinetic energy K = (\u\ ) /2 
and the average square polymer elongation E = (tra). 
The total energy of the system Et = K + ^^E is the 
sum of kinetic and elastic contributions, and its balance 
equation, from (pp)-© is: 



dE 7 
~dt 



(4) 



As it is shown in Fig. [TJ the system does not depart 
from the laminar fixed point up to Wi = 10 which is ap- 
proximately the critical Wi obtained in [18|. For higher 
values of the Weissenberg number, a depletion in the ki- 
netic energy is observed, accompanied by a growth in the 
trace of the conformation tensor, larger than its laminar 
prediction (tra) = 2 + Wi 2 /A. Therefore, polymers drain 
energy from the mean flow in order to elongate and, once 
sufficiently stretched they appreciably start to back-react 
on the flow itself, through the coupling term V • a in 
(|2j). The increase with respect to the laminar growth 
indicates this coupling is very efficient in sustaining the 
stretching of polymers. 
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where ej = (f • u) is the energy input, e v = v((d a up) ) 
is the viscous dissipation, and the last term represents 
elastic dissipation. The balance of forcing and dissipa- 
tion produces a steady state, = in which the total 
energy is constant, at least in a statistical sense (i.e. on 
time scales larger than the typical duration of fluctua- 
tions). 

The mean values of K and E as a function of the Weis- 
senberg number are presented in Fig. [TJ Here and in the 
following the Weissenberg number Wi = Ur/L and the 
Reynolds number Re = UL/v t are defined in terms of 
the amplitude U of the mean velocity, which is measured 
from velocity profiles. Indeed, a remarkable feature of the 
Kolmogorov flow is that even in the turbulent regime the 
mean velocity and conformation tensor are accurately de- 
scribed by sinusoidal profiles [32|, though with different 
amplitudes with respect to the laminar fixed point. 



FIG. 1. (Color online) Mean kinetic energy (K) (a) and mean 
square polymer elongation (E) (b) in statistically stationary 
states as a function of Wi for Re c± 1. Dashed lines are 
the laminar predictions from J3} and error bars represent the 
amplitude of fluctuations. 

We observe here that, due to the present flow geom- 
etry, one of the diagonal elements of the conformation 
tensor is much larger than the other, namely an ^> 022 • 
This implies that the trace of the tensor is almost en- 
tirely given by an, and the same is true for the first 
normal stress difference Ni oc (an — 022) ~ 0"n, so that 
(tra) « (Ni). Together with the growth of the trace of 
the conformation tensor, the growth of the first normal 
stress difference is a clear signature of the elasticity of 
the flow. An inspection at the behavior of the compo- 
nents of the conformation tensor, Fig. [2j indeed reveals 
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that (an) > (022) and, in turn, (cr 22 ) > |(cr 12 )|- In 
particular we find (en) ~ 10(022) independently of Wz, 
and (cr 22 ) « 10|<cr 12 >| for Wi > 10. Although (cr 12 ) and 
((J22) always remain small compared to (<tii), their be- 
havior with Weissenberg number clearly shows a transi- 
tion around Wi ^10, indicating the establishment of a 
nonzero transversal force, growing with Wi, potentially 
capable of sustaining a secondary flow with nonzero com- 
ponent in the direction perpendicular to that of the base 
flow. 




FIG. 2. (Color online) Average values of (011) — ((722), pro- 
portional to the first normal stress difference iVi, and of ((722) 
as a function of Wi. The line corresponds to the Wi 2 behav- 
ior of the laminar fixed point. Inset: ratio (an)/ (022) and 
|(c"i2)| as a function of Wi. 

Let us now discuss the behavior of temporal fluctua- 
tions of the global quantities K (kinetic energy) and E 
(trace of the conformation tensor), when the Weissenberg 
number is increased, which provides further information 
on the change of dynamical regime. In Fig. [3] we report 
K and E versus time, for several values of Wi, above 
the elastic instability threshold Wi c ~ 10. The picture 
emerging from this figure is that of a transition from con- 
stant to randomly fluctuating states, via periodic ones. 

Below Wi = 14, K and E are constant, although at 
different values with respect to the laminar fixed point. 
We remark that the constant behavior does not neces- 
sarily mean the flow is time independent, as it is also 
compatible with the presence of steady, uniformly trav- 
eling patterns of the form A(x — vt) propagating with 
velocity v. 

When Wi > 14 the kinetic energy and the squared 
polymer elongation oscillate around their average values. 
To assess the frequency content of the above fluctuating 
temporal signals, we performed a Fourier analysis, com- 
puting the power spectrum 
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where [0,T max ] is an interval of time in the statistically 



FIG. 3. (Color online) Time behavior of kinetic energy K (a) 
and average square polymer elongation E (b) in statistically 
stationary states, for Re ~ 1 at Wi = 11.1 (continuous red 
lines), Wi = 14 (dashed green lines), Wi = 16.25 (dot-dashed 
blue lines) and Wi = 22 A (dotted pink lines). In (a) Wi 
increases from top to bottom, while in (b) from bottom to 
top, in agreement with the mean values shown in Fig. [T] 



steady state; the results are shown in Fig. |4] (similar re- 
sults are obtained for E). This procedure allowed to 
identify two transitions: at Wi > 14 an explicit time 
dependency appears, signaled by a single frequency peak 
in the spectra of K and E. In this range of parameters 
the system is in a periodic state, characterized by global 
oscillations in time. The peak frequency uu* is found to 
be slightly larger than the frequency uja = 27r /Ta as- 
sociated to the advective time scale, Ta = Lq/U, and 
smaller than that associated to the polymer relaxation 
time uj t = 2tt/t (which is in turn much smaller than the 
frequency ur = 2ir/Tr related to the inverse velocity gra- 
dient Tr = L/U of the mean flow): uja < uo* < uj t ujr- 
The time dependency of the kinetic energy implies that of 
the velocity field u = u(x,t). It is known from the theory 
of dynamical systems that this is a sufficient condition 
for the appearance of Lagrangian chaos in a 2D incom- 
pressible flow and the chaoticity of Lagrangian trajecto- 
ries implies the establishment of mixing features [33| . In 
previous numerical simulations [8| , we measured the La- 
grangian Lyapunov exponent A, defined as the mean rate 
of separation of two infinitesimally close particles trans- 
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ported by the flow. This quantity, which gives an esti- 
mation of the inverse mixing time, was indeed found to 
be positive in a corresponding range of Wi > 10 values, 
implying the onset of chaotic (Lagrangian) dynamics. 
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FIG. 4. (Color online) Frequency spectra Sk(oj) of K(t) at 
Re = 0.875, Wi = 14 (a); Re = 0.8125, Wi = 16.25(b); 
Re = 0.7925, Wi = 19.02 (c). The inset in (c) is the frequency 
spectrum plotted with lin-log axes to show the exponential 
decay (black line). 

The present picture persists at increasing Weissenberg 
numbers, until at Wi = 16.25 a second frequency, asso- 
ciated to a slower time scale, shows up. For Wi > 19.02 
several frequencies are present, indicating more disor- 
dered oscillations. In this regime, the spectral amplitudes 
seem to approximately follow an exponential decay. 

In this latter case, the random fluctuations in time 
of the global quantities point at the onset of irregu- 
lar dynamics associated to a chaotic, mixing, multi- 
scale state. A fixed time snapshot of the vorticity field 
((x,y) = d x u y — d y u x at Wi = 22.4 is shown in Fig. [5] 
where an irregular, disordered pattern is clearly observ- 
able superimposed to the basic periodic flow. 

The statistical characterization of this turbulent-like 
state is given by the spectrum of velocity fluctuations in 
the wavenumber k domain. As shown in Fig. [6j the spec- 
trum displays a power-law behavior E(k) ~ k~ a , with 




FIG. 5. (Color online) Snapshot of the vorticity field C(x, y) = 
V x u(x,y) for Re ~ 0.7 and Wi ~ 22.4. Black (white) 
corresponds to minimum (maximum) vorticity. 



exponent a > 3, reflecting the activity of a whole interval 
of spatial modes. The fast decay of the spectrum, tells 
us that the velocity and velocity gradients are primar- 
ily determined by the integral scale, i.e. on the system 
size. This implies that elastic turbulence corresponds to a 
temporally random, spatially smooth flow, dominated by 
strongly non-linear interactions between few large scale 
modes. 




FIG. 6. (Color online) Spectrum of velocity fluctuations E(k), 
in the wavenumber domain for Re c± 0.7 and Wi ~ 22.4; 
the line corresponds to the best fit k~ 3 ' 8 . The spectrum is 
computed by averaging over several realizations of the velocity 
field. 



B. Wavy patterns 

As discussed in the previous section, for values of the 
Weissenberg number larger than Wi c « 10, we observe 
a transition to non stationary states. We now focus on 
features of such states in space and in time. 

For Wi ~ 11 a typical filamental pattern emerges in 
both velocity (or vorticity) and a fields. Snapshots of 
these fields in the (x, y) plane are shown in Fig. [7] These 
filamental structures are found to perform a wavy motion 
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along the horizontal direction of the mean flow. The pos- 
sibility of observing "elastic waves" in polymer solutions 
was theoretically predicted within a simplified uniaxial 
elastic model [7] which has strong formal analogies with 
MHD equations, but their experimental observation is 
still lacking. 



(a) (b) 




FIG. 7. (Color online) Snapshots of different fields at a given 
time: (a) CO, 2/), (b) <Ju(x,y), (c) a 12 (x,y), (d) cr 22 (x,2/); 
Re — 0.925, Wi = 11.1. In all cases black means minimum 
and white maximum intensity. 

In a narrow range of Wi the patterns are found to 
simply translate without changing their form. Already 
at Wi = 14 temporal oscillations of pattern shapes are 
manifest, particularly on filament tails and in regions be- 
tween them. Nevertheless, wavy motion is still clearly 
detectable. At larger Weissenberg number, filament al 
patterns are significantly less stable. Indeed they are ob- 
served to form and disrupt in the course of time, as they 
travel in both positive and negative horizontal directions. 

Since patterns described so far travel along the longitu- 
dinal direction, it is natural to wonder what their typical 
speed is, how this possibly varies with the transversal 
coordinate, and how these quantities change when the 
elasticity of the solution is increased. We measured the 
speed of these elastic waves by tracking peaks of inten- 
sity in the fields (i,j = 1,2), and following their 
motion along the x— direction (Fig. [8]). 

The results relative to the cases where coherent motion 
of structures was detected are reported in Table U and 
Fig. [9] In all the considered cases we found an average 
wave speed intensity close to half the peak fluid velocity: 
(\v\) ~ U/2. Concerning spatial variations of \v\, these 
show little dispersion S\v\ around the mean at Wi — 11.1, 
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FIG. 8. (Color online) (a) Component an of the conforma- 
tion tensor, (b) Space-time plot obtained by following the 
evolution in time of a horizontal cross section (corresponding 
to the dotted line in (a)) of an(x,y, t). Time is plotted ver- 
sus position; lines of constant slope mean the pattern moves 
at constant velocity. Due to periodic boundary conditions, 
waves leaving the box on the left re-enter on the right. Here 
Re = 0.875, Wi = 14. 



as it can be seen also in the plot (Fig. [9J of the profiles 
along the transversal coordinate y, which are essentially 
constant. The speed of different parts of a wave, for in- 
stance its tip and its tail, tracked by cross-sections at 
different y coordinates, move at approximately the same 
speed, pointing to rigid motion of the structure, which 
travels substantially undeformed. We remark these dy- 
namics are compatible with the constancy of the global 
quantities shown in Fig. [3j due to the fact that K and E 
are integrals over the space domain by definition, and to 
the presence of periodic boundary conditions. When the 
Weissenberg number is increased (Wi = 14 in Fig. [9^l), 
wave speed profiles become less constant, indicating more 
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Wi 


Re 


Wi 


U 


(\v\)±8\v\ 


(\v\)/U 


1 


12 


0.925 


11.1 


3.7 


1.938 ±0.016 


0.52 


1 


16 


0.875 


14.0 


3.5 


1.90 ±0.12 


0.54 


1 


20 


0.813 


16.3 


3.25 


1.70 ±0.27 


0.52 


0.1 


12 


0.093 


11.1 


3.7 


1.704 ±0.022 


0.46 


0.05 


12 


0.049 


11.7 


3.9 


1.722 ±0.032 


0.44 



TABLE I. Measure of the typical speed of waves for various 
values of the Reynolds and Weissenberg numbers. Here Reo 
and Wio are computed from the laminar fixed point velocity 
amplitude Uo =4; U is the measured intensity of the mean 
velocity profile, from which Re and Wi are estimated. 



complex motion. The larger scatter of v values reflects 
the larger uncertainty S\v\ in the measure of the propaga- 
tion speed of patterns, which are now found to swing at 
^/-locations corresponding to their tails, and the presence 
of secondary less steady filamental structures in regions 
between the better defined primary waves. At even larger 
values of the Weissenberg number (Wi = 16.3), waves 
interact considerably more, starting to move in both di- 
rections and to break and reform after they "collide"; 
the distribution of wave speeds now exhibits significantly 
more variability along the y— direction. In fact, at the 
present Reynolds number (Re ~ 1), this is the largest 
Wi at which we succeeded in tracking waves. Above this 
value, patterns form more easily in various parts of the 
spatial domain, but interactions among them dominate 
the dynamics, producing a disordered flow in both space 
and time. 

At the smallest value of Wi — 11.1 we find that the 
dependence of the mean velocity on Re is very weak (see 
Fig. [9)3), confirming the elastic nature of these instabili- 
ties. 

In order to get some insight into the physical mecha- 
nism producing filamental patterns, we now consider the 
spatial structure of the components u x , u y of the velocity 
field, for Re c± 1 and Wi > Wi c . The velocity field can be 
thought as the sum of the base flow and a secondary flow: 
u x = uox + Uy = uo y +u' y \ at the laminar fixed point, 
of course u x = u$ x = Uocos(y/L), u y = uo y = 0. Here 
we are especially concerned with the establishment of a 
nonzero transversal velocity component and the breaking 
of translational invariance along the longitudinal direc- 
tion, above the elastic instability threshold. Snapshots 
at fixed time of the fields u x (x, y, t) and u y (x, y, t), com- 
puted from the vorticity field ( shown in Fig. [T^l, are 
presented in Fig. [TOk fTUb . 

In the horizontal component of velocity u x , the peri- 
odicity of the base flow is clearly visible. The transver- 
sal component u y , on the other hand, displays localized 
quadrupolar patterns of positive and negative values, on 
a background of zero velocity. Just to fix the ideas, con- 
sider, for instance, the pattern centered around y = Lq/2, 
Lq = 2tt being the height of the computational domain, 
on the right part of Fig. [TOa . Notice that the position of 
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FIG. 9. (Color online) Transversal profiles of wave speed: 

(a) measurements at fixed Re c± 1 and varying Weissenberg 
number Wi = 11.1 (red), Wi = 14 (green), Wi = 16.3 (blue); 

(b) measurements at fixed Wi c± 11.1 and varying Reynolds 
number Re = 0.925 (red), Re = 0.093 (green), Re = 0.049 
(blue) . In both plots wave speeds measured in different fields 
(Cj a i2 C22) are presented. Wave speeds are made non- 
dimensional dividing by the intensity of the mean velocity 
profile U, the vertical coordinate is normalized to the box 
size Lq. 




FIG. 10. (Color online) Snapshots in the (x,y) plane of dif- 
ferent fields at Re = 0.925, Wi = 11.1. (a) Secondary trans- 
verse flow u y = Uy] (b) longitudinal flow u x = uo x + u x ; 
(c) square polymer elongation tra. In all cases black means 
minimum and white maximum intensity. Due to the positive- 
definiteness of the conformation tensor, in (c) black corre- 
sponds to zero average square elongation. 



this pattern matches exactly that of the tip of the most 
clearly visible wave in the vorticity field (Fig. [7^). In- 
deed, this structure is found to be stable and to travel 
along the horizontal direction, following the propagation 
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of the vorticity filament. The lobes of these quadrupolar 
patterns (with white, black indicating positive, negative 
values respectively), appear in correspondence to loca- 
tions of strong primary strain, where gradients of the 
base flow attain maximal values in magnitude. In these 
regions polymers are significantly stretched, by the terms 
(Vw) T • a and cr • ( Vu) in eq. © . In Fig. [TUb a snapshot, 
at the same time, of the trace of the conformation ten- 
sor, which measures the square elongation of polymers, 
is reported. Strong elongations appear in thin filamen- 
tal structures localized around y coordinates where u x 
changes sign or, equivalently, where d y u x is maximum. 
As a consequence, it is in this regions that they can exert 
a relevant feedback on the flow, through the elastic stress 
term ^ V • a. 

T 

In other words, the transverse flow develops where 
the feedback of polymers on the velocity dynamics is 
strongest. Finally, the establishment of the transversal 
component of the flow is accompanied by a reduction of 
the longitudinal one, due to incompressibility (Fig. HUb). 
The result of this whole process is the formation of re- 
circulation zones, which are evident in a plot of the cor- 
responding stream function i/j (such that u x = —d y ijj, 
u y = d x ip) (Fig. [TT]) . These in turn result to be orga- 
nized in an array of counter-rotating vortices of equal 
strength, with a basic periodicity of 2 in the longitudinal 
direction, which propagates horizontally. 




X 



FIG. 11. (Color online) Stream function ip(x,y) for Re — 
0.925, Wi = 11.1; black corresponds to minimum, white to 
maximum. 

Concerning the feedback of polymers on the flow field, 
the transversal velocity is found to be sustained by a 
persistent correlation between the transversal compo- 
nents of the velocity field u y and the elastic body force 

fy el ^ = (2r)v/r) (V • <t) , confirming what anticipated in 
the preceding discussion. The ^—profiles of these quanti- 
ties are shown in Fig. [I2j at x ~ 4.7, that is immediately 



on the left of the center of the quadrupolar pattern dis- 
cussed above. Though this figure refers to a single real- 
ization at fixed time, this positive correlation is a feature 
persisting over time. In Fig. [12] it is also shown the pro- 
file of W = fy el)u y/(\fy Cl) \ • Kl), whicn is the 

sign of 

the elastic force's power. Though this fluctuates wildly, 
it is clearly positive in finite size regions located around 
y ~ 7r, corresponding to the extension of the quadrupolar 
pattern on the right in fig [TUh . This observation further 
supports the role of elastic contributions in the establish- 
ment of the transversal flow, and hence in the appearance 
of filamental structures. 
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FIG. 12. (Color online) Transverse polymer body force fy 
and transverse velocity u y (multiplied by a factor 150 for rea- 
sons of visualization) versus the vertical coordinate y. Inset: 
elastic force's power as a function of the vertical coordinate. 
Here Re = 0.925, Wi = 11.1 and the horizontal coordinate is 
x ~ 4.7. This corresponds to a vertical line immediately on 
the left of the center of the quadrupolar pattern in Fig. ITOl 

To conclude, the picture emerging from this descrip- 
tion is that filamental patterns are associated to the for- 
mation of arrays of counter-rotating vortices, sustained 
by the dynamical coupling between velocity and local- 
ized structures where polymers are strongly elongated. 
Indeed this is probably the simplest way to break the 
translational symmetry of the velocity field along the hor- 
izontal direction. The heuristic argument developed here 
is summarized in the scheme in Fig. [T3j As a final re- 
mark, it may be interesting to observe that this allows 
to understand the direction of wave propagation. This 
results from the direction of the base velocity field, from 
which vortices of intensity of order U/2 form. Due to the 
periodicity of the base flow, both positive and negative 
horizontal directions of propagation are possible. The 
selected direction then only depends on where the posi- 
tive correlation between a random perturbation of u y and 

the elastic body force fy Gl ^ happens to be established. As 
shown in Fig. [I4j at increasing the elasticity of the so- 
lution, the probability to find positive correlations gets 
larger, which implies the possibility of forming filaments 
in a larger fraction of the domain. This reflects the ex- 
perimental observation that when Wi is sufficiently large, 
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filaments' propagation may happen in both directions. 



IV. SUMMARY AND CONCLUSIONS 



+ 
+ 




+ 
+ 



FIG. 13. (Color online) Schematic illustration of the forma- 
tion of vortices. Yellow star points indicate locations of max- 
imum velocity gradient in absolute value, where polymers are 
significantly stretched and the transverse flow is more likely 
to occur; small black arrows correspond to the secondary 
transversal flow; big black arrows give the direction of propa- 
gation of waves. The left column contains the sign of uq x and 
the right one that of £. 
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In this work we have studied the dynamics of a two- 
dimensional parallel flow, namely the Kolmogorov shear 
flow, of a viscoelastic fluid, close to the elastic instabil- 
ity threshold predicted in [18] within the framework of 
linear stability analysis. By means of direct numerical 
simulations of the Oldroyd-B model, we investigated the 
destabilization of this flow induced by the elastic forces 
associated to the dynamics of polymer molecules in the 
solution. Above a critical Weissenberg number Wi c « 10 
a transition to new dynamical states was observed. In 
this regime of large elasticity, the system was found to 
be characterized by the presence of filamental structures 
propagating in the direction of the mean flow and global 
oscillations introducing temporal dependencies in the ob- 
servables. This phenomenon allows the establishment of 
mixing features. Close to the critical Weissenberg num- 
ber the structures can be described as traveling waves 
associated to steady patterns of counter-rotating vortices 
sustained by the coupling between elastic forces and ve- 
locity field. The onset of wavy motion of filaments was 
found to occur even in the limit of very small Reynolds 
numbers, typically for Re > 0.05, at Wi > Wi c . 

Increasing the Weissenberg number, the system evolves 
towards more complex dynamics, where filaments are ob- 
served to continuously form and break as a consequence 
of interactions among them. As a consequence, several 
time scales and finer spatial scales result to be involved 
in the dynamical behavior, which displays strongly non- 
linear features resembling those of a turbulent flow. For 
instance, when Wi > 20 (and Re < 1) the power spec- 
trum of velocity fluctuations E(k) is given by a power-law 
k~ a of exponent a larger than 3, in close agreement with 
experimental findings |l|. 

The present results refer to a simplified viscoelastic 
model in an idealized geometrical configuration. Nev- 
ertheless, this simple setup proved useful to detect the 
appearance of coherent structures in the form of "elas- 
tic waves" , driving the transition towards mixing states, 
and to reproduce the basic phenomenology observed in 
experiments. We hope these findings may help in the un- 
derstanding of the phenomenon of elastic turbulence in 
more realistic configurations and in the comprehension 
of the elastic instability mechanism in parallel flows. 



FIG. 14. (Color online) Probability density function of the 
sign of transverse flow's power. At growing Wi positive values 
are more likely to occur. Here the Reynolds number is always 
Re ~ 1. 
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